Water Resources in Copenhagen during 20th century
Before we start: data wrangling
First load your spatial analysis toolkit
library(sf)
library(tidyverse)
library(spatstat)
library(spatialkernel)
library(googlesheets4)
library(leaflet)Next, load your spatial data
Reading layer `bydel' from data source `C:\Users\Adela\Documents\RStudio\1_Teaching\BetweenCultureAndNature\data\bydel.shp' using driver `ESRI Shapefile'
Simple feature collection with 10 features and 4 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 12.45305 ymin: 55.61284 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84
Simple feature collection with 6 features and 4 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 12.46344 ymin: 55.63174 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84
id bydel_nr navn areal_m2
5 6 6 Vanl<f8>se 6699011
6 5 5 Valby 9234177
7 4 4 Vesterbro-Kongens Enghave 8472602
8 1 1 Indre By 10488466
9 9 9 Amager <d8>st 9820410
10 2 2 <d8>sterbro 9858740
geometry
5 MULTIPOLYGON (((12.4982 55....
6 MULTIPOLYGON (((12.52434 55...
7 MULTIPOLYGON (((12.54553 55...
8 MULTIPOLYGON (((12.72897 55...
9 MULTIPOLYGON (((12.63082 55...
10 MULTIPOLYGON (((12.59777 55... [1] 3 7 8 10 6 5 4 1 9 2 [1] "N<f8>rrebro" "Br<f8>nsh<f8>j-Husum"
[3] "Bispebjerg" "Amager Vest"
[5] "Vanl<f8>se" "Valby"
[7] "Vesterbro-Kongens Enghave" "Indre By"
[9] "Amager <d8>st" "<d8>sterbro" Simple feature collection with 10 features and 2 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 12.45305 ymin: 55.61284 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84
navn geometry
1 N<f8>rrebro MULTIPOLYGON (((12.53704 55...
2 Br<f8>nsh<f8>j-Husum MULTIPOLYGON (((12.46894 55...
3 Bispebjerg MULTIPOLYGON (((12.5383 55....
4 Amager Vest MULTIPOLYGON (((12.58271 55...
5 Vanl<f8>se MULTIPOLYGON (((12.4982 55....
6 Valby MULTIPOLYGON (((12.52434 55...
7 Vesterbro-Kongens Enghave MULTIPOLYGON (((12.54553 55...
8 Indre By MULTIPOLYGON (((12.72897 55...
9 Amager <d8>st MULTIPOLYGON (((12.63082 55...
10 <d8>sterbro MULTIPOLYGON (((12.59777 55...
Name
1 N<f8>rrebro
2 Br<f8>nsh<f8>j-Husum
3 Bispebjerg
4 Amager Vest
5 Vanl<f8>se
6 Valby
7 Vesterbro-Kongens Enghave
8 Indre By
9 Amager <d8>st
10 <d8>sterbroNext attach the attribute data
wc <- read_sheet("https://docs.google.com/spreadsheets/d/1iFvycp6M6bF8GBkGjA2Yde2yCIhiy5_slAkGF-RUF7w/edit#gid=0",
col_types = "cnnnnnnnn")
wc# A tibble: 100 x 9
Suburb suburb_id flats wc_access bath year bath_communal_ct wc_communal_ct
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Indre~ 1 16280 11310 3800 1950 NA NA
2 Chris~ 1 5490 3900 900 1950 NA NA
3 Voldk~ 1 13460 12690 4560 1950 NA NA
4 Øster~ 2 30820 28900 13750 1950 NA NA
5 Indre~ 3 28700 24380 5910 1950 NA NA
6 Ydre ~ 3 21710 20410 5800 1950 NA NA
7 Veste~ 4 25850 23930 3730 1950 NA NA
8 Konge~ 4 6270 6240 4240 1950 NA NA
9 Valby 5 14430 13970 8190 1950 NA NA
10 Viger~ 5 7700 7580 5050 1950 NA NA
# ... with 90 more rows, and 1 more variable: hot_water <dbl>Data on washing resources in Copenhagen now looks good and tidy, but its spatial resolution is better than the provided polygons (as in we have multiple rows that all fit within one suburb id). We therefore need to aggregate the data before we attach it to the spatial polygons
wcdata <- wc %>% group_by(year, suburb_id) %>% summarize(flats = sum(flats), bath = sum(bath),
wc_access = sum(wc_access), warmH20 = sum(hot_water), communal_wc = sum(wc_communal_ct),
communal_bath = sum(bath_communal_ct))
wcdata# A tibble: 50 x 8
# Groups: year [5]
year suburb_id flats bath wc_access warmH20 communal_wc communal_bath
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1950 1 35230 9260 27900 8530 NA NA
2 1950 2 30820 13750 28900 10750 NA NA
3 1950 3 50410 11710 44790 8810 NA NA
4 1950 4 32120 7970 30170 6110 NA NA
5 1950 5 22130 13240 21550 10830 NA NA
6 1950 6 10260 6780 10120 6270 NA NA
7 1950 7 27260 14790 26770 13280 NA NA
8 1950 8 19270 15000 18980 14690 NA NA
9 1950 9 23960 12470 22950 11210 NA NA
10 1950 10 18000 9030 16200 7800 NA NA
# ... with 40 more rowsJoin the data with the spatial representations
Now we can join the data with the spatial polygons
Simple feature collection with 50 features and 11 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 12.45305 ymin: 55.61284 xmax: 12.73425 ymax: 55.73271
geographic CRS: WGS 84
First 10 features:
id bydel_nr navn areal_m2 year flats bath wc_access warmH20
1 1 1 Indre By 10488466 1981 26413 14035 22546 NA
2 1 1 Indre By 10488466 1950 35230 9260 27900 8530
3 1 1 Indre By 10488466 1965 32470 12780 32450 NA
4 1 1 Indre By 10488466 1970 30440 11386 22381 NA
5 1 1 Indre By 10488466 1955 33490 9740 33430 9130
6 2 2 <d8>sterbro 9858740 1950 30820 13750 28900 10750
communal_wc communal_bath geometry
1 NA NA MULTIPOLYGON (((12.72897 55...
2 NA NA MULTIPOLYGON (((12.72897 55...
3 9510 2280 MULTIPOLYGON (((12.72897 55...
4 NA NA MULTIPOLYGON (((12.72897 55...
5 NA NA MULTIPOLYGON (((12.72897 55...
6 NA NA MULTIPOLYGON (((12.59777 55...
[ reached 'max' / getOption("max.print") -- omitted 4 rows ] [1] "id" "bydel_nr" "navn" "areal_m2"
[5] "year" "flats" "bath" "wc_access"
[9] "warmH20" "communal_wc" "communal_bath" "geometry" Lets look at the data in a map.
Flats through time
library(tmap)
wc1950 <- wc_spatial %>%
filter(year==1950)
tmap_mode(mode = "view" )
tm_shape(wc_spatial) +
tm_facets(by = "year")+
tm_borders(col = "black",
lwd = 1) +
tm_polygons("flats",
style = "pretty")+
tm_scale_bar(position = c("LEFT", "BOTTOM"),
breaks = c(0, 2, 4),
text.size = 1) +
tm_compass(position = c("RIGHT", "TOP"),
type = "rose",
size = 2) +
tm_credits(position = c("LEFT", "BOTTOM"),
text = "Adela Sobotkova, 2021") +
tm_layout(main.title = "Copenhagen Flats",
#bg.color = "lightblue",
legend.outside = TRUE)library(tmap)
tm_shape(wc_spatial) +
tm_facets(by = "year")+
tm_borders(col = "black",
lwd = 1) +
tm_polygons("flat_per_km",
style = "pretty")+
tm_scale_bar(position = c("LEFT", "BOTTOM"),
breaks = c(0, 2, 4),
text.size = 1) +
tm_compass(position = c("RIGHT", "TOP"),
type = "rose",
size = 2) +
tm_credits(position = c("LEFT", "BOTTOM"),
text = "Adela Sobotkova, 2021") +
tm_layout(main.title = "Copenhagen Flats per sq km",
#bg.color = "lightblue",
legend.outside = TRUE)Access to toilets and baths, per suburb and sq kilometer
Lets calculate the baths and toilets available per square kilometer per each suburb
Continue with communal resources and warm water
Access OSM data for Copenhagen and retrieve ??? (whatever would be relevant??)
The OpenStreetMap contains free and open spatial data for physical features on the ground, with each features’ type being define using key:value pair tags. Each tag describes a geographic attribute of the feature being shown by that specific node, way or relation.
Use:
osmdata:opq()to define the bounding box of the osm requestosmdata:add_osm_feature()to define the key:value pairs you are looking forosmdata:osmdata_sf()to retrieve the osm data.
library(osmdata)
# Create a bounding box
bb <- suburbs %>% st_transform(4326) %>% st_bbox()
plot(bb)
q <- opq(bbox = bb, timeout = 180)
qa <- add_osm_feature(q, key = "amenity", value = "public_bath")
qb <- add_osm_feature(q, key = "amenity", value = "drinking_water")
qc <- add_osm_feature(q, key = "amenity", value = "shower")
qd <- add_osm_feature(q, key = "amenity", value = "toilets")
qe <- add_osm_feature(q, key = "amenity", value = "water_point")
public_bath <- c(osmdata_sf(qa), osmdata_sf(qb), osmdata_sf(qc), osmdata_sf(qd))Clean up OSM data
Use the following code to clean the results and project them in LAEA.
This code:
- removes the duplicated geometries thanks to osmdata::unique_osmdata (see the documentation for details)
- projects into Lambert 93
- keeps the name attribute only
- computes the centroids for the baths stored as polygons
- Eventually, the baths outside CPH suburbs are removed.
bath_uniq <- unique_osmdata(public_bath)
rpoint <- bath_uniq$osm_points %>% filter(!is.na(amenity)) %>% st_transform(32632) %>%
select(name)
rpoly <- bath_uniq$osm_polygons %>% st_transform(32632) %>% select(name) %>% st_centroid()
baths_osm <- rbind(rpoly, rpoint)
baths_osm <- st_intersection(baths_osm, st_transform(suburbs, 32632) %>% st_geometry() %>%
st_union())
# transform also historical baths
baths_cph <- wc_spatial %>% # filter(year==1970) %>%
st_centroid() %>% st_transform(32632) %>% mutate(radius = sqrt(bath_per_km)) %>%
arrange(desc(bath_per_km))Now, let’s display the results in two synchronized mapview maps:
- one with bathing resources in suburbs
- another one with baths extracted from OSM.
- Use the
mapview::syncfunction to display both maps side by side with synchronisation.
library(mapview)
map_osm <- mapview(baths_osm, map.types = "OpenStreetMap", col.regions = "#940000",
label = as.character(suburbs$name), color = "white", legend = FALSE, layer.name = "Baths in OSM",
homebutton = FALSE, lwd = 0.5)
# test map
mapview(baths_cph[, -3], map.types = "Stamen.TonerLite", cex = "radius", legend = FALSE,
col.regions = "#217844", lwd = 0, alpha = 0.4)map_cph <- mapview(baths_cph[, -3], map.types = "OpenStreetMap", col.regions = "#940000",
color = "white", cex = "bath_per_km", legend = TRUE, layer.name = "Baths from 1970",
homebutton = FALSE, lwd = 0.5)
sync(map_osm, map_cph)